﻿using System;
using System.Diagnostics;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Nature.Numerics
{
        /// <summary>
    /// Vector function approximation
    /// </summary>
    public interface IApproximation
    {
        /// <summary>
        /// Finds aproximation of the f-function vector in the vicinity of the x-coordinate
        /// </summary>
        /// <param name="f">function</param>
        /// <param name="x">coordinate</param>
        void Approximate(double[] f, double x);
    }

    /// <summary>
    /// Tabular approximation for the vector function
    /// </summary>
    [Serializable]
    public class Tabulated : IApproximation
    {
        #region Nested Types
        [Serializable]
        struct Point
        {
            /// <summary>
            /// Coordinate
            /// </summary>
            public double X { get; set; }

            /// <summary>
            /// Vector function
            /// </summary>
            public double[] F { get; set; }

            /// <summary>
            /// Vector function derivative defined at the X-coordinate
            /// </summary>
            public double[] DfDx { get; set; }

            /// <summary>
            /// Linear approximation in the vicinity of the X-coordinate
            /// </summary>
            /// <param name="f"></param>
            /// <param name="x"></param>
            public void Approximate(double[] f, double x)
            {
                for (int iva = 0; iva < f.Length; ++iva)
                {
                    f[iva] = this.F[iva] + (x - this.X) * this.DfDx[iva];
                }
            }
        }
        #endregion

        #region Private Fields
        private Point[] m_points;
        private double m_logFactor;
        #endregion

        #region Private Methods

        private static VectorFunctionDelegate ConvertToVectorFunction(IEnumerable<ScalarFunctionDelegate> operatorCollection)
        {
            return
                delegate(double[] f, double x)
                {
                    ScalarFunctionDelegate[] func = operatorCollection.ToArray();
                    for (int iva = 0; iva < func.Length; ++iva)
                    {
                        f[iva] = func[iva](x);
                    }
                };
        }

        private Point FindPoint(double x)
        {
            //x(1) = x(0) * (1 + del)
            //x(2) = x(1) * (1 + del)
            //x(2) = x(0) * (1 + del)**2
            //x(3) = x(0) * (1 + del)**3
            //...
            //x(n) = x(0) * (1 + del)**n
            //(1+del)**n = x(n) / x(0)
            // n * ln(1+del)  = ln(x(n)/x(0))
            // n = ln(x(n)/x(0)) / ln(1+del)            
            double doubleIndex = Math.Log(x / Min) / m_logFactor;
            int ipo = (int)doubleIndex;
            Point p = m_points[ipo];
            if (p.X > x)
            {
                ipo--;
                p = m_points[ipo];
                Debug.Assert(p.X < x, "algorithm failed");
            }
            return p;
        }
        #endregion




        /// <summary>
        /// Number of variables
        /// </summary>
        public int Nva
        {
            get;
            private set;
        }

        /// <summary>
        /// Number of approximation points
        /// </summary>
        public int Npo
        {
            get;
            private set;
        }

        /// <summary>
        /// left hand side boundary of the approximation domain
        /// </summary>
        public double Min
        {
            get;
            private set;
        }


        /// <summary>
        /// right hand side boundary of the approximation domain
        /// </summary>
        public double Max
        {
            get;
            private set;
        }

        /// <summary>
        /// Relative increment
        /// </summary>
        public double Increment
        {
            get;
            private set;
        }

        /// <summary>
        /// Instantinates the tabulated approximation
        /// </summary>
        /// <param name="operatorCollection">collection of scalar functions to be approximated</param>
        /// <param name="min">left hand side boundary of the approximation domain</param>
        /// <param name="max">right hand side boundary of the approximation domain</param>
        /// <param name="delta">relative approximation increment (e.g. 0.01 for the 1% growing)</param>
        public Tabulated(IEnumerable<ScalarFunctionDelegate> operatorCollection, double min, double max, double delta) :
            this(ConvertToVectorFunction(operatorCollection), operatorCollection.Count(), min, max, delta)
        {

        }

        /// <summary>
        /// Instantinates the tabulated approximation
        /// </summary>
        /// <param name="op">vector function</param>
        /// <param name="nva">number of variables</param>
        /// <param name="min">left hand side boundary of the approximation domain</param>
        /// <param name="max">right hand side boundary of the approximation domain</param>
        /// <param name="delta">relative approximation increment (e.g. 0.01 for the 1% growing)</param>
        public Tabulated(VectorFunctionDelegate op, int nva, double min, double max, double delta)
        {
            if (op == null)
            {
                throw new System.ArgumentNullException();
            }
            if (nva < 1)
            {
                throw new System.ArgumentOutOfRangeException();
            }

            Min = min;
            Max = max;
            Nva = nva;
            Increment = delta;
            m_logFactor = Math.Log(1.0 + delta);
            List<double> grid = new List<double>();
            for (double x = min; x <= max; x += delta * x)
            {
                grid.Add(x);
            }
            if (false == grid.Contains(max))
            {
                grid.Add(max);
            }
            Npo = grid.Count;
            m_points = new Point[Npo];
            for (int ipo = 0; ipo < Npo; ++ipo)
            {
                m_points[ipo].X = grid[ipo];
                m_points[ipo].F = new double[Nva];
                op(m_points[ipo].F, m_points[ipo].X);
                m_points[ipo].DfDx = new double[Nva];
            }
            for (int ipo = 0; ipo < Npo - 1; ++ipo)
            {
                Point curr = m_points[ipo];
                Point next = m_points[ipo + 1];
                for (int iva = 0; iva < Nva; ++iva)
                {
                    curr.DfDx[iva] = (next.F[iva] - curr.F[iva]) / (next.X - curr.X);
                }
            }
            m_points.Last().DfDx.Assign(0.0);
        }

        #region IApproximation members
        /// <summary>
        /// Retrieves the linear approximation of the tabulated vector function
        /// </summary>
        /// <param name="f"></param>
        /// <param name="x"></param>
        public void Approximate(double[] f, double x)
        {
            Debug.Assert(f.Length == Nva, "invalid array length");

            if (x < Min)
            {
                f.CopyFrom(m_points[0].F);
                return;
            }
            if (x > Max)
            {
                f.CopyFrom(m_points.Last().F);
                return;
            }

            Point p = FindPoint(x); ;

            for (int iva = 0; iva < Nva; ++iva)
            {
                f[iva] = p.F[iva] + (x - p.X) * p.DfDx[iva];
            }

        }
        #endregion
    }

}
